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1.0 Introduction 


The program PLANECNG calculates velocities for Earth-to-Moon or Moon-to-Earth trajecto¬ 
ries. The flight to be analyzed originates in a circular orbit of any inclination and altitude about 
one of the bodies, and culminates in a circular orbit of any inclination and altitude about the 
other body. An intermediate AV and plane change occurs at the Lunar sphere of influence (SOI), 
the region where the vehicle is near its lowest velocity in the trajectory, and therefore where it is 
able to make the plane change with the lowest AV. 

A given flight may penetrate the SOI at a number of points. Each point has associated with it a 
unique set of AV’s and total velocity. This program displays the velocities, in matrix form, for a 
representative set of SOI penetration points. An SOI point is identified by projecting Lunar 
latitude and longitude onto the SOI. The points reported for a given flight are defined by the 
user, who provides a starting longitude and latitude, and an increment for each. A matrix is built 
with ten longitudes forming the columns and 19 latitudes forming the rows. This matrix is 
presented in six different reports, each report containing different velocity or node information in 
the body of the matrix. 

Section 2.0 of this document describes the inputs provided by the user to define the flight profile. 
Section 3.0 describes the contents of the six reports that are produced as outputs. Section 4.0 
contains the instructions to execute this program. The appendices provide a look at the structure 
and details of the program code. 
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2.0 Program Inputs 

The following paragraphs discuss the inputs provided by the user. The prompt is the message 
displayed by the program onto the screen. The input variable is the variable assigned to the 
user’s response. The description provides information about how to respond to the prompt. 


1. Prompt: INPUT OUTBOUND OR RETURN 

Input variable: MD 

Description: Enter "OUTBOUND" for Earth-to-Moon trajectory calculations. Enter 
"RETURN" for Moon-to-Earth trajectory calculations. Be sure to type the 
double quotation marks. 

2. Prompt: INPUT NODE OPTION 1 OR 2 

Input variable: NP 

Description: Given that the Earth-Moon trajectory will intercept the SOI at a particular 
latitude and' longitude, and given a desired destination circular orbit 
inclination (greater than or equal to the SOI latitude), there exist two 
destination orbits that provide a solution (see figure 2-1). Node option I 
directs the program to concern itself with the orbit in which the SOI 
longitude is less than 90° away from the destination orbit’s ascending 
node. Node option 2 directs the program to concern itself with the orbit in 
which the SOI longitude is greater than 90° away from the destination 
orbit’s ascending node. 

3. Prompt: INPUT PERIGEE ALTITUDE OF EARTH ORBIT (NMI) 

Input variable: HPE 

Description: Enter the height above the Earth’s surface of the Earth circular orbit, in 
nautical miles. 

4. Prompt: INPUT PERIGEE ALTITUDE OF LUNAR ORBIT (NMI) 

Input variable: HPM 

Description: Enter the height above the Lunar surface of the Lunar circular orbit, in 
nautical miles. 

5. Prompt: INPUT EARTH DEPARTURE JULIAN DATE 

Input variable: TIMJ 

Description: Enter the origin trans-SOI injection date (Earth departure date for 
outbound trajectories, and Lunar departure date for return trajectories) in 
Julian day format, where January 1, 2000 is day 2,451,545. Refer to 
Section C of "The Astronomical Almanac of the Year 1988". Day 
2451545 is the default value if zero is entered in this field. 

6. Prompt: INITIAL LONGITUDE 

Input variable: ALONI 

Description: Enter initial sphere of influence longitude for the output matrices. This 
value will become the heading for column 1 of the matrices. 
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Figure 2-1 




Prompt: INPUT INCREMENT FOR THE MAP 

Input variable: DELLON 

Description: Enter longitude increments for the output matrices. Applied to the initial 
longitude, this value defines the subsequent column headings of the 
matrices. Longitudes for outbound trajectories should be between zero 
and -90 degrees. Longitudes for return trajectories should be between 
zero and +90 degrees. 


8. Prompt: INPUT INITIAL LATITUDE 

Input variable: ALATI 

Description: Enter initial sphere of influence latitude for the output matrices. This 
value will become the heading for row 1 of the matrices. 

9. Prompt: INPUT INCREMENT FOR THE MAP 

Input variable: DELLAT 

Description: Enter latitude increments for the output matrices. Applied to the initial 
latitude, this value defines the subsequent row headings of the matrices. 

10. Prompt: INPUT EARTH ORBIT TO LUNAR ORBIT INCLINATION 

Input variable: AINCEO 

Description: Enter Earth circular orbit inclination, in degrees. This is not what is 
typically considered inclination (i.e., a measurement taken from the 
Earth’s equatorial plane), but rather the angle between the plane of the low 
Earth orbit and the plane of the Moon’s orbit about the Earth. 

11. Prompt: INPUT LUNAR ORBIT TO LUNAR ORBIT INCLINATION 

Input variable: AINCM 

Description: Enter Lunar circular orbit inclination, in degrees. This is the angle 

between the plane of the low Lunar orbit and the plane of the Moon’s orbit 
about the Earth. 

12. Prompt: INPUT FLIGHT TIME 

Input variable: FTIM 

Description: Enter the desired total flight time from origin trans-SOI injection to 
destination circular orbit injection, in hours. 
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3.0 Program Outputs 


This section describes the contents and shows a sample of the six reports generated by the 
program. 


Report #1 

Velocity Map For Inbound/Outbound Trajectories 

The top section of this report repeats the input values entered by the user. The second section is 
a 10 X 19 matrix of total velocities required to fly the profile described by the inputs. Each cell 
corresponds to a particular latitude and longitude on the Lunar sphere of influence. A plane 
change bum may occur at any one of these coordinates, and the value of the corresponding cell is 
the total velocity required for the flight if the bum occurs at that location. It is the sum of the 
transfer orbit injection AV, sphere of influence plane change AV, and destination circular orbit 
injection AV. The axes of the matrix are determined by the initial longitude, latitude, and 
increments specified by the user. 

The third section of the report is a summary of key data corresponding to the matrix cell 
containing the lowest total velocity. This data includes: 

• X-, Y-, and Z-components of velocity at the sphere of influence just before and just 
after the plane change bum (VX, VY, VZ). 

• Total magnitude of the velocity at the sphere of influence just before and just after the 
plane change bum (VEL). 

• Flight path angle at the sphere of influence just before and just after the plane change 
bum (GAMA). 

• Azimuth of the sphere of influence point from the Earth and from the Moon (AZM). 

• Earth and Lunar orbit inclinations (AINC). 

• Transfer trajectory and circular orbit insertion ascending/descending nodes (ANODE). 

• Earth-to-SOI and SOI-to-Moon times of flight (TIME). 

• Transfer trajectory and circular orbit insertion AV’s (DVCIR). 

■ Sphere of influence AV (DVPHER). 

• Total AV (DVTOTAL). 
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VELOCITY MAP FOR OUTBOUND TRAJECTORIES NODE OPTION, 2 
DATE 7-SEP-88 TIME 08:35:41 

JUILIAN DAY 2451545. PERIGEE ALT (NMI) EARTH = 250. MOON = 60. 

TRANSLUNAR FLIGHT TIME (HR) = 60.0 INCL EARTH 30.0 INCL MOON = 30. 
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Report #2 

Delta Velocity At Sphere Of Influence Map For Translunar/Transearth Trajectories 

This report is a matrix of delta velocities that occur at the sphere of influence. Each cell 
corresponds to a particular latitude and longitude on the sphere of influence at which a plane 
change bum may occur. 
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DELTA VELOCITY AT SPHERE OF INFLUENCE MAP FOR TRANSLUNAR TRAJECTORIES NODE OPTION 

DATE 7-SEP-88 TIME 08:35:41 

TRANSLUNAR FLIGHT TIME (HR) = 60.0 INCL EARTH 30.0 INCL MOON = 30.0 
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Report #3 

Ascending Node Of Lunar Orbit Map For Translunar/Transearth Trajectories 

This report is a matrix of ascending nodes at the Moon. Each matrix cell displays an ascending 
node corresponding to a trajectory that penetrates the SOI at the stated SOI latitude and 
longitude. The node describes the Lunar longitude at which the SOI-Moon trajectory passes 
through the Earth-Moon plane during Lunar circular orbit insertion for outbound flights, or 
during trans-SOI injection for return flights. 
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Report #4 

Ascending Node Of Earth-Moon Transfer Trajectory Map For Translunar/Transearth Trajecto¬ 
ries 

This report is a matrix of ascending nodes at the Earth. Each matrix cell displays an ascending 
node corresponding to a trajectory that penetrates the SOI at the stated SOI latitude and 
longitude. The node describes the Earth longitude at which the Earth-SOI trajectory passes 
through the Earth-Moon plane during Earth circular orbit insertion for return flights, or during 
trans-SOI injection for outbound flights. 
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ASCENDING NODE OF EARTH-MOON TRANSFER TRAJECTORY MAP FOR TRANSLUNAR TRAJECTORIES 
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Report #5 

Delta Velocity At Earth For Earth-Moon Transfer Trajectory Map For Translunar/Transearth 
Trajectories 

This report is a matrix of trajectory velocities at the Earth. Each cell corresponds to a particular 
targeted latitude and longitude on the SOI at which a plane-change bum may occur. For Earth- 
to-Moon trajectories, the values in these cells are the trans-SOI insertion trajectory AV’s. For 
Moon-to-Earth trajectories, the values in these cells are the Earth circular orbit insertion AV’s. 
This program assumes a departure from and arrival to circular orbit. 
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DELTA VELOCITY AT EARTH FOR EARTH-MOON TRANSFER TRAJECTORY MAP FOR TRANSLUNAR TRAJECTORIES 

NODE OPTION 2 

DATE 7-SEP-88 TIME 08:35:41 

TRANSLUNAR FLIGHT TIME (HR) = 60.0 INCL EARTH 30.0 INCL MOON = 30.0 
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Report #6 

Delta Velocity At Moon For Earth-Moon Transfer Trajectory Map For Translunar/Transearth 
Trajectories 


This report is a matrix of trajectory velocities at the Moon. Each cell corresponds to a particular 
targeted latitude and longitude on the SOI at which a plane-change bum may occur. For Earth- 
to-Moon trajectories, the values in these cells are Lunar circular orbit insertion AV’s. For Moon- 
to-Earth trajectories, the values in these cells are the trans-SOI insertion trajectory AV‘s. 

Lunar circular orbits at the specified altitudes are assumed. 
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DELTA VELOCITY AT MOON FOR EARTH-MOON TRANSFER TRAJECTORY MAP FOR TRANSLUNAR TRAJECTORIES 

NODE OPTION 2 

DATE 7-SEP-88 TIME 08:35:41 

TRANSLUNAR FLIGHT TIME (HR) = 60.0 INCL EARTH 30.0 INCL MOON = 30.0 


o 

o 

o 

o 

♦ 

o 

o 

o 

r* 

r* 

r- 

r* 

r- 

r- 

o 

o 

o 

o 

♦ 

o 

o 

t 

o 

• 







(O 

CO 

CO 

CO 

CO 

CO 

o 







o 







VO 

VO 

VO 

VO 

VO 

VO 

VO 







CD 

1 







CM 

CM 

CM 

CM 

CM 

CM 

CM 







o 

o 

# 

o 

o 

o 

o 

o 

CO 

CM 

CM 

CM 

CM 

CO 

in 

o 

o 

o 

o 

o 

o 

• 







r- 

rH 

H 

rH 

rH 

r- 

CO 







o 







VO 

r- 

r- 

r- 

r- 

VO 

VO 







00 

1 







CM 

CM 

CM 

CM 

CM 

CM 

CM 







o 

o 

♦ 

o 

o 

o 

o 

o 

o 

rH 

o 

o 

o 

CD 

o 

o 

o 

o 

o 

o 

o 

. 







m 

CD 

CD 

CD 

<D 


r- 







o 







r- 

r- 

r- 

r- 

r- 

r* 

VO 







r- 







CM 

CM 

CM 

CM 

CM 

CM 

CM 







o 

o 

o 

o 

o 

o 

o 



00 

00 

CO 

o 

00 

o 

o 

o 

o 

o 

o 

. 







r- 

r- 

rH 

rH 

r- 

CO 

VO 







o 







00 

00 

CD 

CD 

00 

00 

VO 







VO 

1 







CM 

CM 

CM 

CM 

CM 

CM 

CM 







o 

o 

o 

o 

o 

o 

o 

CM 

rH 

r- 

r- 

rH 

rH 

r- 

o 

o 

o 

o 

o 

o 

. 







r- 

VO 

o 

o 

VO 

r- 

00 







o 







00 

CD 

o 

o 

CD 

00 

r- 







to 

1 







CM 

CM 

CO 

CO 

CM 

CM 

CM 







o 

o 

o 

o 

o 

o 

o 

in 

m 

o 

CD 

CM 


in 

o 

o 

o 

o 

o 

o 

. 







rH 

o 

o 


in 

rH 

00 







o 







CD 

o 

rH 

rH 

o 

CD 

r- 














CM 

CO 

CO 

CO 

CO 

CM 

CM 







o 

o 

o 

o 

o 

o 

o 

CO 


r* 

00 

00 

CM 

CO 

o 

o 

o 

o 

o 

o 

• 







rH 

o 



<D 

rH 

00 







o 







<d 

o 

rH 

CM 

o 

CD 

r- 







CO 

1 







CM 

CO 

CO 

CO 

CO 

CM 

CM 







o 

o 

o 

o 

o 

o 

o 

00 

VO 

CD 

r** 

CM 

r- 

rH 

o 

o 

o 

o 

o 

o 

• 







VO 

m 


CD 

o 

vo 








o 







00 

CD 

o 

o 

o 

00 

r- 







04 

1 







CM 

CM 

CO 

CO 

CO 

CM 

CM 







o 

o 

o 

o 

o 

o 

o 

rH 

VO 

CD 

CD 

m 

o 

o 

o 

o 

o 

o 

o 

o 

• 







00 

CO 

o 

o 

VO 

00 

o 







o 







r* 

CO 

CD 

CD 

00 

r* 

r- 







rH 

1 







CM 

CM 

CM 

CM 

CM 

CM 

CM 







o 

o 

o 

o 

o 

o 

o 

CD 

00 

00 

r- 

00 

CD 

CM 

o 

o 

o 

o 

o 

o 

• 







CD 

CO 

CO 

CO 

CO 

in 

CM 







o 







VO 

r- 

r* 

r- 

r- 

VO 

VO 














CM 

CM 

CM 

CM 

CM 

CM 

CM 








A 


o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

£ 


• 

• 

• 

• 

* 

* 

• 

• 

• 

♦ 

* 

• 

• 

• 

• 

• 

* 

• 

• 

o 


o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

£ 

g 







00 

CM 

rH 


rH 

1 

CM 

1 

CO 

1 








16 














4.0 Program Execution Instructions 


The following instructions describe the steps to be taken by the user to execute this program. 
A. Obtain access to the DEC VAX minicomputer and sign on with user identification. 


B. 

C. 

D. 

E. 


F. 


At the $ prompt, RTTN PLANECNG. 

When prompted by the program, enter the program inputs. See section 2.0 for a 
discussion of the inputs. 


After the last input has been entered, the program will execute for approximately 10 to 30 
minutes. When it is completed, the message "FORTRAN STOP will appear. 


The program outputs will be placed in a file named PLANECNG. OUT;### where### is 
a system generated version number of the report. To print the reports, type the following 

at the $ prompt: 

tvpf pi ANECNG.OUT;### 


To re-execute the program with new parameters, begin again at step (B) above. 
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Appendix A - Program Flow Chart 
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Declare the program 
matrices 

















Record the current 
date and time 



Initialize current a V 
hold variables 
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Increment column 
counter by 1 


Increment longitude 








Have x 
all 19 matrix 
<ows been processe 

\ ? / 









24 
























Identify the angle 
between the Earth-to-Moon 
line and Earth-toSOI line 


Determine the Earth-based 
latitude and longitude 
of the SOI 


Determine the Earth 
ascending node for the orbit 
In which the trans-SOI bum 
will occur 



Identify the angle between 
the Moon-to-Earth line 
and the Moorvto-SOI line 


Determine the azimuth of 
the SOI point from 
the Moon 
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GAMACALC 
(determine time-of- 
fllght #1) 



GAMACALC 
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Calculate A V 
at Lunar perigee 



POSVELMO 

















State minimum total A V 
and corresponding SOI 
longitude, latitude, and 
velocity 



Store minimum SOI A V 
and corresponding SOI 
longitude, latitude, and 
velocity 
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Set Indicators to show 
hyperbolic orbit. Print 
screen message announcing 
hyperbolic orbit. 










S 



Calculate circular 
velocity at perigee 














B 


Calculate true anomaly 


Calculate eccentric anomaly 


Is orbital radius 
sat SOI > apogee radius^ 


Display orbital/apogee radius 
difference, SOI lat. and long., 
Q, semi-major axis, eccen¬ 
tricity, flight path angle, and 
perigee velocity 










POSVELMO 
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Determine the radius of 
the Lunar orbit from the 
Earth-Moon baricenter 


Calculate the Moon's 
velocity in the direction 
of Its orbit 






MOON 



3 ' 





























Appendix B - Program Code Listing 
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u u u 


c 

c 

c 


*** JACK FUNK’S PLANE CHANGE PROGRAM *** TO 87-57 TASK 1.4 
MAIN 


IMPLICIT REAL* 16 (A-H,0-Z) 

CHARACTER* 10 MD, HEAD, TRAJ, TRAJM, TRAJE 


CHARACTER* 8 TEMP 
CHARACTER*9 DATP 

DIMENSION DELV(10,19), VELMOUT(10,19),ALONO(10),ALATP(19), 

* PAGEl(lO), PAGE2(10,19), PAGE3(10,19) ,PAGE4(10,19), 

* PHIM(2), PAGE5(10,19), PAGE6(10,19) 


OPEN OUTPUT FILE 


C 


OPEN (UNIT = 1 .FILE = ’PLANECNG.OUT’ .STATUS = ’NEW’) 
OUTPUT TO FILE; IP: OUTPUT TO SCREEN; IS: OUTPUT TO DEBUG 


IP = 1 
IS =5 

DPR =57.29578 
PI =3.1415926535 
CMUE = 1.407647E+16 
CMUM = 1.731432E+14 
FTNM =6076.115 
REE =20925741. 

REMO = 5.7039E+6 
RREM = 207559 * FTNM 
FTM =.3048 

PRINT *, ’INPUT OUTBOUND OR RETURN ’ 
READ *, MD 

IF (MD .EQ. ’RETURN’ )THEN 
MOOD = -1 

HEAD = ’TRANSEARTH’ 


ELSE 
MOOD =1 

HEAD = ’TRANS LUNAR’ 

ENDIF 

PRINT *,’INPUT NODE OPTION 1 OR 2 ’ 

READ *, NP 

0 PRINT *, ’INPUT PERIGEE ALTITUDE OF EARTH ORBIT (NMI) ’ 

READ * HPE 

PRINT *, ’INPUT PERIGEE ALTITUDE OF LUNAR ORBIT (NMI) ’ 

READ *, HPM 

RPE = HPE * FTNM + REE 

RPM = HPM * FTNM + REMO 

PRINT *, ’INPUT EARTH DEPARTURE JULIAN DATE ’ 

READ *’, TEMJ 

IF (QABS (TEMJ ) .LE. 0.001 ) TIMJ = 2451545. 

PRINT *, ’LONGITUDE FOR OUTBOUND TRAJECTORIES SHOULD BE 
* BETWEEN 0 AND -90 DEG’ 


40 



PRINT *, ’AND RETURN TRAJECTORIES BETWEEN 0 AND +90 DEG’ 
PRINT *, ’INITIAL LONGITUDE’ 

READ *, ALONI 

PRINT *, ’INPUT INCREMENT FOR MAP ’ 

READ *, DELLON 

PRINT *, ’INPUT INITIAL LATITUDE’ 

READ *, ALATI 

PRINT *, ’INPUT INCREMENT FOR MAP ’ 

READ *, DELL AT 

PRINT *, ’INPUT EARTH ORBIT TO LUNAR ORBIT INCLINATION ’ 
READ *, AINCEO 
AINCEO = AINCEO / DPR 

PRINT *, ’INPUT LUNAR ORBIT TO LUNAR ORBIT INCLINATION ’ 
READ *, AINCM 
AINCM = AINCM /DPR 
C PAGE-1 

AINC = AINCM 

AINCE = AINCEO 

PRINT *, ’INPUT FLIGHT TIME ’ 

READ *, FTIM 
CALL DATE(DATP) 

CALL TIME(TIMP) 

15 CONTINUE 

WRITE (IP,7) MD, NP 

7 FORMAT (T12,’ VELOCITY MAP FOR ’,A8,’ TRAJECTORIES 
* NODE OPTION Ml) 

WRITE (IP,17) DATP, TIMP 
17 FORMAT (T27,’DATE’,A10,’ TIME’,A10) 

WRITE (IP,27) TTMJ, HPE, HPM 

27 FORMAT (/,’ JUILIAN DAY ’, F8.0,’ PERIGEE ALT (NMI) EARTH = ’ 
*F4.0, ’ MOON = ’, F4.0) 

WRITE (IP,37) FTIM, AINCEO * DPR, AINCM * DPR 
37 FORMAT (’ TRANS LUNAR FLIGHT TIME (HR) = ’,F5.1 ’ INCL 
♦EARTH ’ F5.1, ’ INCL MOON = ’, F5.1) 

IPRINT = 0 
II =1 
NN = 1 
DVMIN =99999. 

YDLO = QSQRT (CMUE / RREM) 

ALAT = ALATI/DPR 
ALON = ALONI/DPR 
VELM = VELMI 
DVSIM = 99999. 

21 CONTINUE 
VEL = VELM 
DELV(NN,H) =99999. 

PAGEl(NN) =0.0 
C PAGE-2 

PAGE2(NN,II) = 0.0 



PAGE3(NN,II) = 0.0 
PAGE4(NN,II) = 0.0 
PAGE5(NN,II) = 0.0 
PAGE6(NN,II) = 0.0 
C PRINT*,” 

C PRINT*,’ ALAT ALON VEL DVT DVSI MOON EARTH 
C *TIMEE TIMEM RREM’ 

22 CONTINUE 
C 

C CALCULATE DATA FOR POINT IF LATITUDE IS LESS THAN INCLINATION 
C LOAD ZEROS IN DELTA VEL 

IF(QABS(QSIN(ALAT )) .LE. QABS (QSIN(AINCM))+ .00001) THEN 
DO 8 NV = 0, 40 

VELM = 800. + 100. * QFLOAT(NV ) 

VEL = VELM 
ICALL = 1 
GOTO 25 

1000 CONTINUE 

IF (QABS(TIEM) .LE. .0000000000000001) GOTO 23 
C STORE VALUES FOR OUTPUT 

IF(DVTOTAL .LT. DELV(NN,II)) THEN 
DELV (NNJI) = DVTOTAL 
PAGE1 (NN) = DVTOTAL 
PAGE2 (NNJI) = DEL VEL 
VELMOUT (NN,II) = VELM 
C PAGE-3 

ALATP (II) = ALAT * DPR 
PAGE3 (NNJI) = ANODEM * DPR 
PAGE4 (NNJU) = ANODEE * DPR 
PAGE5 (NN,II) = DVCIRE 
PAGE6 (NN,II) = DVCIRM 

C WRITE (IS.57) ALAT * DPR, ALON * DPR, VELMOUT (NN,E) 

C * ,DELV(NN,II), PAGE2(NN,II), TRAJM, TRAJE, TIEM, 

C * TIMM, RREMER 

C 57 FORMAT (1X,F6.1,1X,F6.1,2X,F7.0,1X,F8.1,1X,F7.1,1X, 

C * A5,1X,A5,1X,F6.1,1X,F6.1,1X,F8.3) 

C WRITE (IS,67) VXE, VXM, VYE, VYM, VZE, VZM 

C 67 FORMAT (’ VXE ’,F7.1,’ VXM ’,F7.1,’ VYE \F7.1,’ VYM ’ 

C * F7.1, ’ VZE ’,F7.1, ’VZM \F7.1) 

ENDIF 

8 CONTINUE 
ELSE 

PAGE1 (NN) =0.0 
DELV (NN,13) = 0.0 
PAGE2 (NNJI) = 0.0 
PAGE3 (NNJI) = 0.0 
PAGE4 (NNJI) = 0.0 
PAGE5 (NNJI) = 0.0 
PAGE6 (NNJI) = 0.0 
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ENDIF 

23 CONTINUE 

ALONO (NN) = ALON * DPR 

IF (NN .EQ. 10 .AND. IPRINT .EQ. 0 ) THEN 

WRITE (IP,77) ALONO(l), ALONO(2), ALONO(3), ALONO(4), 

* ALONO(5), ALONO(6), ALONO(7),ALONO(8),ALONO(9),ALONO(10) 

77 FORMAT (/,’ ALON> \ 10(2X,F5.1 )) 

WRITE (IP,87) 

87 FORMAT (’ALAT’) 

ENDIF 

IF (NN .LT. 10) THEN 

ALON = (ALONI + DELLON * QFLOAT (NN )) / DPR 
NN = NN + 1 
GOTO 21 
ENDIF 
IPRINT = 1 

WRITE (IP,97) ALAT * DPR, PAGE1(1),PAGE1(2), 
*PAGE1(3),PAGE1(4),PAGE1(5),PAGE1(6), 
*PAGE1(7),PAGE1(8),PAGE1(9),PAGE1(10) 

97 FORMAT (IX, F5.1.2X, 10(1X,F6.0)) 

IF (II .LT. 19) THEN 

ALAT = (ALATI + DELLAT * QFLOAT (II)) / DPR 

n =11+1 

NN =1 

ALON = ALONI / DPR 
GOTO 21 
ENDIF 

VELM = VELMMIN 
ALAT = ALATMIN 
ALON = ALONMIN 
ICALL = 2 
GOTO 25 
2000 CONTINUE 

WRITE (IP,107) NP , T1 , 

107 FORMAT (///,T23,’MINIMUM VELOCITY PRINT NODE OPTION ,11) 

WRITE (IP,117) DATP, TIMP 
117 FORMAT (T27,’DATE’A10,’TIME’A10) 

WRITE (IP, 127) FTTM, ALAT * DPR, ALON * DPR, AINCEO * DPR, 

* AINCM * DPR 

127 FORMAT (/,’ FLIGHT TIME = ’,F4.0,’ LAT = ’,F6.1,’ LON = ’,F6.1, 

* ’ LUNAR AINC= ’,F5.1,’ EARTH AINC = ’,F5.1) 

WRITE (IP, 137) 

137 FORMAT (’BODY VX VY VZ VEL GAMA AZM 
♦AINC ANODE TIME DVCIR’) 

WRITE (IP, 147) VXM,VYM,VZM,VMAGM,GAMAM * DPR, 

♦AZMM ♦ DPR, AINCM ♦ DPR, ANODEM * DPR, TIMM, DVCIRM, DELVEL 
147 FORMAT (’ MOON ’,4(1X,F6.0),1X,F5.1,1X,F6.1,2X,F5.1,2X,F6.1, 

* 1X,F5.1,1X,F7.0,’ DVPHER = ’,F7.0) 

WRITE (IP,1475) VXE, VYE, VZE, VELE, GAMAE * DPR, AZME DPR, 
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* AINCE * DPR, ANODEE * DPR, TIEM, DVCIRE, DVTOTAL 
1475 FORMAT (’ EARTH’,4( 1X,F6.0), 1X,F5.1,1 X,F6.1,2X,F5.1,2X,F6.1, 

* 1X,F5.1,1X,F7.0,’ DVTOTAL = ’, F7.0) 

WRITE (IP.157) CHAR(12) 

157 FORMAT (’ ’AO 

WRITE (IP, 167) HEAD,NP 

167 FORMAT (T10,’DELTA VELOCITY AT SPHERE OF INFLUENCE MAP 
*FOR ’,A10,’ TRAJECTORIES NODE OPTION ’ Jl) 

WRITE (IP,177) DATP.TIMP 
177 FORMAT (T27,’DATE’A10,’ TIME’AIO) 

WRITE (IP, 187) HEAD,FTIM, AINCEO * DPR, A1NCM * DPR 
187 FORMAT (’ ’,A10,’ FLIGHT TIME (HR) = \F6.1,’ INCL EARTH ’, 

*F6.1,’ INCL MOON =’,F5.1) 

WRITE (IP, 197) ALONO(l), ALONO(2), ALONO(3), ALONO(4), 
*ALONO(5),ALONO(6), ALONO(7), ALONO(8), ALONO(9), ALONO(IO) 
197 FORMAT (/,’ ALON> ’,10(2X, F5.1)) 

18 CONTINUE 
WRITE (IP,207) 

207 FORMAT (’ALAT’) 

DO 28 NPI = 1 , 19 

WRITE (IP,217) ALATP(NPI), PAGE2(1,NPI),PAGE2(2,NPI), 

* PAGE2(3,NPI),PAGE2(4,NPI)JPAGE2(5,NPI),PAGE2(6J S JPI), 

* PAGE2(7,NPI),PAGE2(8,NPI),PAGE2(9,NPI),PAGE2(10,NPI) 

217 FORMAT(1X,F5.1,1X,10(1X,F6.1)) 

28 CONTINUE 

WRITE (IP, 227 ) CHAR(12) 

227 FORMAT (’ ’, Al) 

WRITE (IP, 237 ) HEAD , NP 

237 FORMAT (T6,’ASCENDING NODE OF LUNAR ORBIT MAP FOR ’,A10, 

* ’ TRAJECTORIES NODE OPTION ’,11) 

WRITE (IP, 247 ) DATP, TIMP 

247 FORMAT (T27,’DATE’A10,’ TIME’,A10) 

WRITE (IP, 257 ) HEAD, FTIM, AINCEO * DPR, AINCM * DPR 
257 FORMAT (’ ’,A10,’ FLIGHT TIME (HR) = ’,F6.1,’ INCL EARTH’, 
*F6.1,’ INCL MOON =’,F6.1) 

WRITE (IP,267) ALONO(l), AJLONO(2), ALONO(3), ALONO(4), 
*ALONO(5),ALONO(6), ALONO(7), ALONO(8), ALONO(9), ALONO(IO) 
267 FORMAT (/,’ ALON> ’, 10(2X, F5.1)) 

WRITE (IP,277) 

277 FORMAT (’ ALAT’) 

DO 48 NPI = 1 , 19 

WRITE (IP,287) ALATP(NPI),PAGE3( 1 ,NPI),PAGE3(2,NPI), 

* PAGE3(3,NPI),PAGE3(4,NPI),PAGE3(5,NPI),PAGE3(6JSfPI), 

* PAGE3(7,NPI), PAGE3(8,NPI),PAGE3(9,NPI),PAGE3(10,NPI) 

287 FORMAT (’ ’,F5.1,1X,10(2X,F5.1)) 

48 CONTINUE 

WRITE (IP, 297) CHAR(.12) 

297 FORMAT (’ ’,A1) 

WRITE (IP, 307 ) HEAD,NP 
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307 FORMAT (T6,’ ASCENDING NODE OF EARTH-MOON TRANSFER TRAJECTORY 
♦MAP FOR \A10,’ TRAJECTORIES NODE OPTION ’,11) 

WRITE (IP, 317 ) DATP.TIMP 
317 FORMAT (T27,’DATE’,A10,’ TIME’,A10) 

WRITE (IP, 327 ) HEAD, FTIM, AINCEO * DPR, AINCM * DPR 
327 FORMAT (IX, A10,’ FLIGHT TIME (HR) = ’,F5.1,’ INCL EARTH ’, 

* F5.1,’ INCL MOON = \F5.1) 

WRITE (IP, 337) ALONO(l),ALONO(2),ALONO(3),ALONO(4),ALONO(5), 
*ALONO(6), ALONO(7), ALONO(8), ALONO(9), ALONO(IO) 

337 FORMAT (/,’ ALON> ’, 10(2X, F5.1)) 

WRITE (IP,347) 

347 FORMAT (’ ALAT’) 

DO 58 NPI = 1 , 19 

WRITE (IP,357) ALATP(NPI),PAGE4(1,NPI),PAGE4(2,NPI), 

* PAGE4(3,NPI),PAGE4(4,NPI),PAGE4(5,NPI),PAGE4(6,NPI), 

* PAGE4(7,NPI),PAGE4(8,NPI),PAGE4(9,NPI),PAGE4(10,NPI) 

357 FORMAT(1X,F5.1,2X,10(1X,F6.2)) 

58 CONTINUE 

WRTTE(IP,367) CHAR(12) 

367 FORMAT(’ ’,A1) 

WRTTE(IP,377) HEAD, NP 

377 FORMAT(T2,’DELTA VELOCITY AT EARTH FOR EARTH-MOON TRANSFER 

* TRAJECTORY MAP FOR ’,A10,’ TRAJECTORIES NODE OPTION ’,11) 
WRTTE(IP,387) DATP, TIMP 

387 FORMAT(T27,’DATE ’, A10,’ TIME ’, A10) 

WRTTE(IP,397) HEAD, FTIM, AINCEO * DPR, AINCM * DPR 
397 FORMAT(’ \A10,’ FLIGHT TIME (HR) = \F6.1,’ INCL EARTH’, 

*F6.1,’ INCL MOON = \F6.1) 

WRTTE(IP,407) ALONO(l), AJLONO(2), ALONO(3), ALONO(4), 

*ALONO(5), ALONO(6), ALONO(7), ALONO(8), ALONO(9), ALONO(IO) 

407 FORMAT!/,’ ALON> ’, 10(2X, F5.1)) 

WRTTE(IP,417) 

417 FORMAT!’ ALAT’) 

DO 68 NPI = 1, 19 

WRITE(IP,427) ALATP(NPI), PAGE5(1,NPI), PAGE5(2,NPI), 

* PAGE5(3,NPI), PAGE5(4,NPI), PAGE5(5,NPI), PAGE5(6,NPI), 

* PAGE5(7,NPI), PAGE5(8,NPI), PAGE5(9,NPI), PAGE5(10,NPI) 

427 FORMAT (1X,F5.1,2X,10(1X,F6.0)) 

68 CONTINUE 

WRITE(IP,437) CHAR(12) 

437 FORMAT(’ ’,A1) 

WRTTE(IP,447) HEAD, NP 

447 FORMAT(T2,’DELTA VELOCITY AT MOON FOR EARTH-MOON TRANSFER 

* TRAJECTORY MAP FOR \A10,’ TRAJECTORIES NODE OPTION ’,11) 

WRITE(IP,457) DATP, TIMP 

457 FORMAT(T27,’DATE’,A10,’ TIME’,A10) 

WRITEUP,467) HEAD, FTIM, AINCEO * DPR, AINCM * DPR 
467 FORMAT(1X,A10,’ FLIGHT TIME (HR) =’,F5.1,’ INCL EARTH’, 

*F5.1,’ INCL MOON =’,F5.1) 
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WRITE(IP,477) ALONO(l), ALONO(2), ALONO(3), ALONO(4), ALONO(5) 

* ALONO(6), ALONO(7), ALONO(8), ALONO(9), ALONO(IO) 

477 FORMAT^’ ALON> 10(2X, F5.1)) 

WRrTE(IP,487) 

487 FORMATS ALAT’) 

DO 78 NPI = 1, 19 

WRITE(IP,497) ALATP(NPI), PAGE6(1,NPI), PAGE6(2,NPI), 

* PAGE6(3,NPI), PAGE6(4,NPI), PAGE6(5,NPI), PAGE6(6,NPI), 

* PAGE6(7,NPI), PAGE6(8,NPI), PAGE6(9,NPI), PAGE6(10,NPI) 

497 FORMAT( 1 X,F5.1,2X, 10(1 X,F6.0)) 

78 CONTINUE 
GOTO 3000 
25 CONTINUE 


COSALAT = QCOS (ALAT ) 

SINALAT = QSIN (ALAT ) 

COSALON = QCOS (ALON ) 

SINALON = QSIN (ALON ) 

COSPHI = COSALAT * COSALON 
COS AINC = QCOS (AINCM / DPR ) 

SINAINC = QSIN (AINCM / DPR ) 

RMR = COSPHI +QSQRT (COSPHI ** 2. - (1. - CMUE / CMUM )) 

RRM = RREM / RMR 

XXM = RRM * COSALAT * COSALON 

XX = -XXM + RREM 

YY = -RRM * COSALAT * SINALON 

ZZ = RRM* SINALAT 

YYM = -YY 


C 

C 


WRITE (IS ,507) XX,YY,ZZ 

507 FORMATS XYZ POSITION AT SPHERE ’,F11.0,1X FI 1 0 IX FI 1 
RRE = QSQRT(XX**2.+YY**2.+ZZ**2) ' ’ 

COSANG A = XX / RRE 


0) 


SINANGA = QSQRT (YY ** 2. + ZZ ** 2.) / RRE 
ANGA = QATAN (SINANGA / COSANG A ) 

ALONX = QATAN (YY / XX ) 

ALATX = QATAN (ZZ / QSQRT (XX ** 2. + YY **2 ) ) 
ANODEE = ALONX - QASIN (QTAN (-ALATX ) / QTAN (AINCI 
END SPHERE OF INFLUENCE CALC 


AAMIN = (RPM + RRM ) / 2. 

EEMIN = (RRM - RPM ) / (RPM + RRM ) 

WMIN = QSQRT (CMUM *(1. -EEMIN ) / RRM ) 

C WRITE(IS,517) WMIN 
C 517 FORMAT (’ MINIMUM VEL AT SI = \F8.4) 
AAMINE = (RPE + RRE ) / 2. 

EEMINE = (RRE - RPE ) / (RPE +RRE ) 

VVMINE = QSQRT (CMUE * (1. - EEMINE ) / RRE ) 
30 CONTINUE 


)) 


CALL GAMACALC(RPM,VELM,RRM,CMUM,COSGAM,VPM,VCIRMP TIMM 
*TRAJM,DPR,ALAT,ALON,FTNM) 

GAMAM = -QFLOAT (MOOD ) * QATAN (QSQRT (l.-COSGAM**2.)/ COSGAM) 
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TEMP = ZZ/RRM/QSIN(AINCM) 

IF (ABS(TEMP) .GT. 1.) THEN 
IF (TEMP .GT. 1.) THEN 
PHIM(l) = PI/2 
ELSE 

PHIM(l) = -PI/2 
ENDIF 
ELSE 

PHIM(l) = QASIN (ZZ/RRM/QSIN(AINCM)) 

ENDIF 

IF (ZZ EQ. 0.) THEN 
PHIM(2) = PI 
ELSE 

PHIM(2) = (PI - PHIM (1 ) * ZZ/QABS(ZZ)) * ZZ/QABS(ZZ) 

ENDIF 

ANODEM = ALON-AVM 

ip/ANODEM .LT.O.) ANODEM = ANODEM + 2.*PI 

CALL VELTRANS(VELM,GAMAM,AZMM,ALAT,ALON, VXM, VYM.VZM, VMAGM,- 

DPR) 

VXM = VXM+XDLO 
VYM = VYM+YDLO 
QQEMIN = 2.1 *RPE/(RRE+RPE) 

VELEM1N = QSQRT(QQEMIN*CMUE/RRE) 

IF(VELE .LT. VELEMIN) VELE = VELEMIN 

TIEMS = TIEM . . . 

C TIMEES CHANGED TO TIEMS, TIMEE CHANGED TO TIEM 

35 CONTINUE 

CALL GAMACALC(RPE,VELE,RRE,CMUE,COSGAME,VPE,VCIRE, 

♦TIME l.TRAJE,DPR, ALAT,ALON,FTNM) 

VELE2 — VELE +10. 

CALL GAMACALC(RPE,VELE2JRRE,CMUE,COSGAME,VPE,VCIRE, 

*TIME2,TRAJE,DPR,ALAT ,ALON,FTNM) 

DELVELE = 10./(TIME2-TIME1)*(FTIM-TIME1-TIMM) 

IF(DELVELE .LT. 500.) THEN 
VELE =VELE + DELVELE 
ELSE 

VELE = VELE + DEL VELE/QABS (DEL VELE)*500. 

CALL F GAMACALC(RPE,VELE,RRE,CMUE,COSGAME,VPE,VCIRE,TIEM, 

♦TRAJE,DPR,ALAT,ALON,FTNM) 

IF(TIEM .EQ. 0.) GOTO 60 
TIMT = TIEM+TIMM 

GAMAE^QFL^AT(MOOD)*CATAN(QSQRT(l.-COSGAME**2.)/COSGAME) 

ALONE = 180. / DPR + ALONX 

ALATE = QATAN(ZZ/QSQRT(XX**2.+YY**2.)) 

50 CONTINUE 
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AINCE = AINCEO 

C ’CALCULATION OF EARTH ORBIT AZIMUTH AT SPHERE OF INFLUENCE 
IF(AINCE .NE. 0.0) THEN 
PHIE = PI-QASIN(ZZ/(RRE*SIN(AINCE)» 

AZME = QATAN2(1./QTAN(A1NCE),QC0S(PHIE)) 

ELSE 

AZME = PI/2 
END IF 
C 

CALL VELTRANS(VELE,GAMAE,AZME,ALATE,ALONE,VXE,VYE,VZE,VELE, 
* DPR) 

DELVEL = QSQRT((VXE-VXM)**2.+(VYE-VYM)**2.+(VZE-VZM)**2.) 
DVCIRM = VPM - VCERMP 
DVCIRE = VPE - VCIRE 

DVE = QSQRT (VCIRE ** 2. + VPE**2.-2. * VCIRE*VPE* 

♦COS (AINCEO - AINCE )) 

DVTSAV = DVTOTAL 
DVTOTAL = DELVEL + DVCIRM + DVCIRE 
C CALCULATE POSITION AND VELOCITY OF MOON 
C ITTERATE FOR FLIGHT TIM TO SPHERE OF INFLUENCE 
TIM = (TIMJ-2451545. )/36525 .+TIEM/876600. 

CALL POSVELMO(TIM,RREMER,XDLO,YDLO) 

RREM = RREMER*20295741. 

IF(ABS(TIEM-TIEMS) .GT. .1) GOTO 25 
IF(DVTOTAL.LT. DVMIN)THEN 
DVM1N = DVTOTAL 
ALONMIN = ALON 
ALATMIN = ALAT 
VELMMIN=VELM 
ENDIF 

IF(DELVEL .LT. DVSIM)THEN 
DVSIM = DELVEL 
ALONSIM = ALON 
ALATSIM = ALAT 
VELMSIM = VELM 
ENDIF 

60 CONTINUE 

IF (ICALL.EQ.l) GOTO 1000 
IF (ICALL.EQ.2) GOTO 2000 
3000 STOP 
END 


SUBROUTINE XYZPOS(RRXALAT,ALON,XX,YX,ZX) 
IMPLICIT REAL * 16 (A-Z) 

XX = -RRX*COS(ALAT)*COS(ALON)+RREM 
YX = -RRX*COS(ALAT)*SIN(ALON) 

ZX= RRX*SIN(ALAT) 

RETURN 


48 



uuuuuuuuuu 


END 


SUBROUTINE POSVELMO(TIM,RRM,XDLO,YDLO) 
IMPLICIT REAL * 16 (A-H,0-Z) 

DELT = 0.5/36525 ./24./3600. 

T1 = TIM-DELT 

CALL MOON(T 1 ,RAM,DECM,RM 1) 

T2 = TIM+DELT 

CALL MOON(T2,RAM,DECM,RM2) 

XDLO = (RM2-RM1)*20925741. 

RRM = (RM2+RMl)/2. 

RRM = RRM 

RRMB = RRM -7.412789E-01 
YDLO = 200570.2/RRMB 
RETURN 
END 


SUBROUTINE MOON (T, RAM, DECM, RM) 

FINDS LOCATION OF MOON IN EQUATORIAL COORDS. AT ANY TIM 
REF: ’87 ASTRONOMICAL ALMANAC 

T IS JULIAN CENTURIES SINCE YEAR 2000 
LAM IS MOON’S ECLIPTIC LONGITUDE 
BETA IS MOON’S ECLIPTIC LATITUDE 
PIE IS HORIZONTAL PARALLAX 
RM IS DIST. TO MOON IN EARTH RADII 
RAM IS RT. ASCENSION OF MOON 
DECM IS MOON’S DECLINATION 
SD IS SEMIDIAMETER OF MOON’S ORBIT 

IMPLICIT REAL * 16 (A-Z) 

C PRINT *, ’ MOON’ 

P = 3.1415926535 
C = P /180 

LAM = C*218.32+C*481267.883*T+C* 6.29 * QSIN(C * 134.9 + C * 
*477198.85 * T) - C * 1.27 * QSIN(C * 259.2 - C * 413335.38 * 

*T) + c * .66 * QSIN(C * 235.7 + c * 890534.23*T) 

LAM = LAM + c * .21 * QSIN(c * 269.9 + c * 954397.7*T) - c * 

*.19 * QSIN(c * 357.5 + c * 35999.05 * T) - c * .11 * 

*SIN(c * 186.6 + c * 966404.05*T) 

beta = c*5.13*QSIN(c*93.3 + c * 483202.03 * T) + c * 

* .28 * QSIN(C * 228.2 + C * 960400.87 * T) -c*.28*QSIN 
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*(c*318.3+c* 60003.18*T)-c*.17*QSIN(c*217.6-c*407332.2 * T) 

pie = c*.9508+c*.0518*COS(c* 134.9 + c * 477198.85 *T) + 
*c*.0095*QCOS(c*259.2-c*413335.38*T)+c* .0078*COS(c * 23 
*5.7+c*890534.23*T)+c*.0028*QCOS(c*269.9+c*954397.7* T) 

SD = .2725 * pie 
RM = 1. /QSIN(pie) 

1 = QCOS(beta) * QCOS(LAM) 

M = .9175 * QCOS(beta) * QSIN(LAM) - .3978 * QSIN(beta) 
n = .3978 * QCOS(beta) * QSIN(LAM) + .9175 * QSIN(beta) 

RAM = QATAN2(M, 1) 

DECM = QASIN(n) 

RETURN 

END 

SUBROUTINE GAMACALC(RPX,W,RRX,CMUX,COSGAMX,VPX,VCIRX, 

*timx,traj,dpr,alatalon,ftnm) 

IMPLICIT REAL * 16 (A-H, O-Z) 

CHARACTER* 10 TRAJ 
IHYPER = 1 
TRAJ = ’ELIPT’ 

QQX = RRX* W**2/CMUX 
IF(QQX-2 .LT. 1.0E-06) QQX = QQX - 1.0E-06 
IF (QQX .GT. 2.) THEN 
IHYPER = -1 
TRAJ = ’HYPER’ 

C PRINT *,’ ***** TRAJECTORY IS HYPERBOLIC ’ 

ENDIF 

AAX = RRX/(2.-QQX) 

IF (AAX. GT . 1.0E12 .OR. AAX. LT. -1.0E12) AAX = -1.0E12 

EEX = l.-RPX/AAX 

PPX = AAX*(1. -EEX ** 2.) 

COSGAMX =QSQRT(RPX/RRX*(1+EEX)/QQX) 

GAMAX = QACOS(COSGAMX) 

VPX = QSQRT(CMUX*(1+EEX)/RPX) 

VCIRX = QSQRT(CMUX/RPX) 

COSTHETAX = (PPX/RRX-Q/EEX 
THETAX = QACOS(COSTHETAX) 

COSAEX = (EEX+COSTHET AX)/( 1+EEX*COSTHETAX) 

ERRRP = RRX-AAX* (1+EEX) 

IF(ERRRP .GT. 0. .AND. AAX .GT. 0.) THEN 
WRITE (IS,537) ERRRP/6076.1155, ALAT*DPR, ALON*DPR 
537 FORMAT (’ RADIUS > APOGEE BY NMI ’,F7.5, F8.0, F7.5) 

WRITE (IS,547) QQX, AAX/FNTM, EEX, GAMAX*DPR, VPX 
547 FORMAT (’ QQX = ’,F7.5,’ AAX = ’,F8.0,’ EEX = \F7.5, 

* ’ GAMAX =’,F5.1,’VPX =’,F7.1) 

ENDIF 
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IF(QQX .GT. 2.) GOTO 101 
C CALC FLIGHT TIM FOR ELIPTICAL ORBITS 
IF (ERRRP .GT. 0.) THEN 
TIMX = 0. 

GOTO 199 
ENDIF 

AEX = QACOS(COSAEX) 

SINAEX = QSQRT( 1-COS AEX**2.) 

TIMX = QSQRT(AAX ** 3. / CMUX)*(AEX-EEX*SINAEX)/3600. 
GOTO 199 
101 CONTINUE 

C CALC FLIGHT TIM FOR HYPERBOLIC ORBITS 
COSHF = COS AEX 
SINHF = QSQRT(COSHF**2.-l.) 

FFX = QLOG(COSHF+QSQRT(COSHF**2.-l.)) 

TIMX = QSQRT(-1. *AAX*AAX*AAX/CMUX)*(EEX*SINHF-FFX)/3600. 
199 CONTINUE 
RETURN 
END 


SUBROUTINE VELTRANS(VEL,GAMA, AZM,ALAT,ALON,VXX,VYX,VZX,VMAG, 
* DPR) 

IMPLICIT REAL * 16 (A-Z) 

RRD = VEL * QSIN (GAMA ) 

RPHID = VEL * QCOS (GAMA) 

VLON = RPHID * QSIN (AZM ) 

VLAT = RPHID * QCOS (AZM ) 

VXR = -RRD * QCOS (ALAT ) * QCOS (ALON ) 

VYR = -RRD * QCOS (ALAT ) * QSIN (ALON) 

VZR = RRD * QSIN (ALAT ) 

VXLA = VLAT * QSIN (ALAT ) * QCOS (ALON ) 

VYLA = VLAT * QSIN (ALAT ) * QSIN (ALON ) 

VZLA = VLAT * QCOS (ALAT) 

VXLO = VLON * QSIN (ALON ) 

VYLO = -VLON * QCOS (ALON ) 

VZLO = 0.0 

VXX = VXR + VXLA + VXLO 
VYX = VYR + VYLA + VYLO 
VZX = VZR + VZLA + VZLO 

VMAG = QSQRT (VXX **2. + VYX ** 2.+ VZX ** 2. ) 

RETURN 

END 
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Appendix C - Program Variables 
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INPUT 

VARIABLE DESCRIPTION 


AINCEO Earth orbit inclination (degrees). This is the angle between the plane of the low 
Earth orbit and the plane of the Moon’s orbit about the Earth. 

AINCM Lunar orbit inclination (degrees). This is the angle between the plane of the low 
Lunar orbit and the plane of the Moon’s orbit about the Earth. 

ALATI Initial Sphere of Influence longitude for map (degrees) 

ALONI Initial Sphere of Influence latitude for map (degrees) 

DELLAT Incremental latitude for map (degrees) 

DELLON Incremental longitude for map (degrees) 

FIIM Flight time for trajectory (hours) 

HPE Holding perigee altitude of Earth orbit (nautical miles) 

HPM Holding perigee altitude of Lunar orbit (nautical miles) 

MD Leg of trip for which to perform calculations (outbound or return) 

NP Node option. Selection of "1" directs the program to concern itself with destina¬ 

tion circular orbits in which the SOI longitude is less than 90° away from the 
destination orbit’s ascending node. Selection of "2" directs the program to 
concern itself with destination circular orbits in which the SOI longitude is 
greater than 90° away from the the destination orbit’s ascending node. 

TIMJ Earth departure Julian date (where January 1, 2000 is day 2,451,545. Refer to 

Section C of "The Astronomical Almanac of the Year 1988"). 

CONSTANT VALUE DESCRIPTION _ 


c 

n/i 80 

Degrees per radian (deg./rad.) 

CMUE 

1.407647E+16 

Gravitational parameter of the Earth (fiVsec 2 ) 

CMUM 

1.731432E+14 

Gravitational parameter of the Moon (ft 3 /sec 2 ) 

DPR 

57.29578 

Degrees per radian (deg./rad.) 

FTM 

0.3048 

Meters per foot (m/ft) 
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FTNM 

6,076.115 

P 

3.1415927 

PI 

3.141593 

REE 

20,925,741 

REMO 

5,703,900 

RREM 

1,261,152,353 


Feet per nautical mile (fit/nmi) 

II (dimensionless) 

II (dimensionless) 

Radius of the Earth (ft) 

Radius of the Moon (ft) 
Earth-Moon distance of centers (ft) 
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VARIABLE 

AAMIN 

AAMINE 

AAX 

AEX 

AINC 

AINCE 

ALAT 

ALATE 

ALATMIN 

ALATP(-) 

ALATSIM 

ALATX 

ALON 

ALONE 

ALONMIN 

ALONO(-) 

ALONSIM 

ALONX 

ANGA 

ANODEE 

ANODEM 

AVM 


DESCRIPTION _ 

Semi-major axis of the SOI-to-Moon trajectory 

Semi-major axis of the Earth-to-SOI trajectory 

Semi-major axis of one of the transfer orbits 

Eccentric anomaly of one of the transfer orbits 

Lunar orbit inclination 

Earth orbit inclination 

Initial latitude for map (radians) 

Latitude of SOI point, measured from the Earth 
Latitude of SOI point associated with the minimum total AV 
ALAT in degrees 

Latitude of SOI point associated with the minimum SOI AV 
Latitude of the SOI point, measured from the Earth 
Initial longitude for map (radians) 

Longitude of SOI point, measured from zero longitude at Earth (-X direction) 
Longitude of SOI point associated with the minimum total AV 
Variable that contains the print matrix column headings (longitudes) 

Longitude of SOI point associated with the minimum SOI AV 
Longitude of the SOI point, measured from the Earth-Moon line at the Earth 
The angle between the Earth-to-Moon line and the Earth-to-SOI-point line. 
Longitude of the Earth ascending or descending node 
Longitude of the Lunar ascending or descending node 

Angle between the Lunar node and the projection of the SOI point onto the Earth- 
Moon plane (v) 


55 




AZM 

AZME 

AZMM 

BETA 

CMUX 

COSAEX 


Azimuth of the SOI point 
Azimuth (from the Earth) of the SOI point 
Azimuth (from the Moon) of the SOI point 
Moon’s ecliptic latitude 
Earth or Moon gravitational parameter 
Cosine of the eccentric anomaly of one of the transfer orbits 
COSAINC Cosine of the Lunar orbit inclination 
COSALAT Cosine of the SOI point latitude 
COSALON Cosine of the SOI point longitude 

COSANGA Cosine of the angle between the Earth-to-Moon line and the Eaith-to-SOI-point 
line 

COSGAM Cosine of the flight path angle at SOI of the SOI-to-Moon trajectory 
COSGAME Cosine of the flight path angle at SOI of the Earth-to-SOI trajectory 
COSGAMX Cosine of the flight path angle at SOI of one of the transfer orbits 

Hyperbolic cosine of the eccentric anomaly of one of the transfer orbits 
Cosine of the angle between the Moon-to-Earth line and the Moon-to-SOI-point 


COSHF 

COSPHI 


line 


COSTHETAXCosine of the true anomaly of one of the transfer orbits 


DATP 

DECM 

DELT 

DELV(v) 

DELVEL 

DELVELE 


Today’s date 
Declination of the Moon 
A fraction of TIM that represents a half-second 
Total AV for pairs of SOI longitudes and latitudes 
AV at the SOI point 

Estimated AV required for the Earth-to-SOI trajectory in order to satisfy the time 
of-flight requirement 
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DVCIRE 

DVCIRM 

DVE 

DVMIN 

DVSIM 

DVTOTAL 

DVTSAV 

EEMIN 

EEMINE 

EEX 

ERRRP 

FFX 

GAMA 

GAMAE 

GAMAM 

GAMAX 

HEAD 

ICALL 

n 

IHYPER 

IP 

IPRINT 

IS 

L 


AV between Earth circular orbit and Earth-to-SOI trajectory 

AV between Lunar circular orbit and Moon-to-SOI trajectory 

Unused plane change variable 

Hold variable for minimum total AV 

Hold variable for minimum SOI AV 

Total AV for flight 

Temporary storage for total AV 

Eccentricity of SOI-to-Moon trajectory 

Eccentricity of Earth-to-SOI trajectory 

Eccentricity of one of the transfer orbits 

Difference between orbital range at SOI and apogee range 

Hyperbolic eccentric anomaly 

Flight path at SOI of one of the transfer orbits 

Flight path at SOI of Earth-to SOI trajectory 

Flight path at SOI of SOI-to-Moon trajectory 

Flight path angle at SOI of one of the transfer orbits 

Heading that indicates leg of journey ("TRANSEARTH" or "TRANSLUNAR") 
Flag to track occurance of call to in-program subroutine. 

Counter (1-19), representing increments in latitude for the output matrices 

Indicator that describes whether an orbit is hyperbolic 

Output number indicating output to file 

Flag for print block that contains matrix column headers 

Output number indicating output to screen 

A geocentric direction cosine 
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LAM 

M 

MOOD 

N 

NN 

NPI 

NV 

PAGE1(») 

PAGE2(v) 

PAGE3(v) 

PAGE4( V ) 

PAGE5(v) 

PAGE6(v) 

PHIE 

PHIM(») 

PIE 

PPX 

QQEMIN 

QQX 

RAM 


Moon’s ecliptic longitude 

A geocentric direction cosine 

Leg of journey (-1: return. +1: outbound) 

A geocentric direction cosine 

Counter (1-10) representing increments in longitude for the output matrices 

Iteration counter. Used to iterate through 19 matrix rows during report printing 

Iteration counter. Used to iterate through 41 values of VELM for each lat¬ 
itude/longitude combination, in search of the minimum AV. 

Total AV for a given SOI longitude and latitude. Cell values for Report #1 
matrix. 

SOI AV for a given SOI longitude and latitude. Cell values for Report #2 matrix. 

Longitude for the Lunar ascending or descending node for a given SOI longitude 
and latitude. Cell values for Report #3 matrix. 

Longitude for the Earth ascending or descending node for a given SOI longitude 
and latitude. Cell values for Report #4 matrix. 

AV between Earth circular orbit velocity and trans-Lunar trajectoiy velocity for a 
given SOI longitude and latitude. Cell values for Report #5 matrix. 

AV between Lunar circular orbit velocity and trans-Earth trajectory velocity for a 
given SOI longitude and latitude. Cell values for Report #6 matrix. 

ANGA + 180°. Used to determine the azimuth of the SOI point where 0° Earth 
longitude is in the -X direction. 

Angle between the Moon-to-Earth line and the Moon-to-SOI-point line 
Horizontal parallax 

Semi-latus rectum of one of the transfer orbits 
Vis-viva parameter for the Earth-to-SOI-point trajectory 
Vis-viva parameter for an orbit 
Right ascension of the Moon 
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RM 

RMR 

RM1 

RM2 

RPE 

RPHID 

RPM 

RPX 

RRD 

RRE 

RREMER 

RRM 

RRMB 

RRX 

SD 

SINAEX 

SINAINC 

SINALAT 

SINALON 

SINANGA 

SINHF 

T 

TEMP 


Distance from the Earth to the Moon in Earth radii 
Ratio of Earth-Moon distance to Moon-SOI distance 
Distance from the Earth to the Moon in Earth radii 
Distance from the Earth to the Moon in Earth radii 
Distance from Earth’s center to Earth perigee orbit 
Orbital path component of the velocity vector 
Distance form Moon’s center to Lunar perigee orbit 
Distance from Earth or Moon center to perigee 
Radial component of the velocity vector 
Distance from the Earth to the SOI point 
Moon’s distance from the Earth (Earth radii) 

Distance from the Moon to the SOI point, when the Moon is at its mean distance 
from the Earth 

Distance of the Moon from the Earth-Moon baricenter 
Distance from Earth or Moon to SOI 
Semi-diameter of the Moon’s orbit 

Sine of the eccentric anomaly of one of the transfer orbits 
Sine of the Lunar orbit inclination 
Sine of the latitude of the SOI point 
Sine of the longitude of the SOI point 

Sine of the angle between the Earth-to-Moon line and the Earth-to-SOI point line 
Hyperbolic sine of the eccentric anomaly of one of the transfer orbits 
Number of Julian centuries since the year 2000 AD 
Temporary hold variable 
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THETAX 

TIEM 

TIEMS 

TIM 

TIM El 

TIME2 

TIMM 

TIMP 

TIMT 

TTMX 

TRAJ 

TRAJE 

TRAJM 

T1 

T2 

VCIRE 

VCIRMP 

VCIRX 

VEL 

VELE 

VELE2 

VELEMIN 

VELM 


True anomaly of one of the transfer orbits 
Earth-to-SOI time of flight (seconds) 

Temporary storage for Earth-to-SOI time of flight 

Time of arrival at SOI point from Earth, in centuries since the year 2000 AD 
Earth-to-SOI time of flight (seconds) 

Earth-to-SOI time of flight (seconds) 

Time of flight from SOI to Lunar perigee 

Time Now 

Total time of flight 

Time of flight from Earth or Moon perigee to SOI 

Text that describes whether an orbit is hyperbolic or elliptical 

Text that describes whether the Earth-to-SOI trajectory is hyperbolic or elliptical 

Text that describes whether the SOI-to-Moon trajectory is hyperbolic or elliptical 

One-half second before TIM 

One-half second after TIM 

Velocity of Earth circular orbit 

Velocity of Lunar circular orbit 

Velocity of the Earth or Lunar circular orbit 

Velocity at SOI of one of the transfer orbits 

Velocity at SOI of the Earth-to-SOI trajectory 

Ten feet per second more than VELE 

Minimum velocity required such that apogee of the trajectory is just at SOI 
Velocity at SOI of the SOI-to-Moon trajectory 
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VELMMIN Velocity at SOI of the SOI-to-Moon trajectory associated with the minimum total 
AV 


VELMOUT(v) Velocity at SOI of the SOI-to-Moon trajectory, for a given SOI longitude and 
latitude 


VELMSIM 

VLAT 

VLON 

VMAG 

VMAGM 

VPE 

VPM 

VPX 

W 

WMLN 

WMINE 

VXE 

VXLA 

VXLO 

VXM 

VXR 

VXX 

VYE 

VYLA 

VYLO 


Velocity at SOI of the SOI-to-Moon trajectory associated with the minimum SOI 
AV 

Latitude component of the velocity vector 

Longitude component of the velocity vector 

Velocity vector magnitude 

Velocity at SOI of the SOI-to-Moon trajectory 

Perigee velocity of the Earth-to-SOI trajectory 

Perigee velocity of the SOI-to-Moon trajectory 

Perigee velocity for one of the transfer orbits 

Velocity at SOI of one of the transfer orbits 

Velocity at SOI point (apogee) of SOI-to-Moon trajectory 

Velocity at SOI point (apogee) of Earth-to-SOI trajectory 

X-coordinate of velocity vector at SOI for Earth-to-SOI trajectory 

X-component of the latitude component of the velocity vector 

X-component of the longitude component of the velocity vector 

X-coordinate of velocity vector at SOI for SOI-to-Moon trajectory 

X-component of the radial component of the velocity vector 

Total X-component of the velocity vector 

Y-coordinate of velocity vector at SOI for Earth-to-SOI trajectory 
Y-component of the latitude component of the velocity vector 
Y-component of the longitude component of the velocity vector 
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VYM Y-coordinate of velocity vector at SOI for SOI-to-Moon trajectory 

VYR Y-component of the radial component of the velocity vector 

VYX Total Y-component of the velocity vector 

VZE Z-coordinate of velocity vector at SOI for Earth-to-SOI trajectory 

VZLA Z-component of the latitude component of the velocity vector 

VZLO Z-component of the longitude component of the velocity vector 

VZM Z-coordinate of velocity vector at SOI for SOI-to-Moon trajectory 

VZR Z-component of the radial component of the velocity vector 

VZX Total Z-component of the velocity vector 

XDLO X-coordinate of the velocity of the Moon 

XX Distance in the X-direction from the Earth to the SOI point’s X-coordinate 

XXM Distance in the X-direction from the Moon to the SOI point’s X-coordinate 

YDLO Y-coordinate of the velocity of the Moon 

YY Distance in the Y-direction from the Earth-Moon line to the SOI point’s Y- 

coordinate 

ZZ Distance in the Z-direction from the Earth-Moon plane to the SOI point 
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Appendix D - Detailed Program Description 
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This section is a step-by-step description of the lines of code in program PLANECNG. Refer to 
Appendix B for the actual program listing, and figures J-l through J-3 for supporting illustra¬ 
tions. 


1. Declare the matrices DELV, VELMOUT, PAGE1, PAGE2, PAGE3, PAGE4, PAGE5, 
PAGE6, PHIM, ALONO, and ALATP. 

2. Open the output files. 

3. Define the program constants. 

4. Read the program inputs. 

a. MD (leg of trip for which to perform calculations; "outbound" or "return") 

b. NP (destination node option) 

c. HPE (perigee altitude of Earth orbit) 

d. HPM (perigee altitude of Lunar orbit) 

e. TIMJ (Earth departure Julian date. Default is 2451545, representing Jan. 1, 2000 
AD) 

f. ALONI, DELLON (Initial longitude and increment for output map) 

g. ALATI, DELLAT (Initial latitude and increment for output map) 

h. AINCEO (Angle between the plane of the low Earth orbit and the Earth-Moon 
plane) 

i. AINCM (Angle between the plane of the low Lunar orbit and the Earth-Moon 
plane) 

j. FILM (flight time) 

5. Calculate the distance from the Earth’s center to Earth perigee orbit. 

RPE = HPE * FTNM + REE 

6. Calculate the distance from the Moon’s center to Lunar perigee orbit. 

RPM = HPM * FTNM + REMO 

7. Store the inclinations in temporary variables (AENCE and AINC). 

8. Record the current date and time (DATP and LIMP). 
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9. 


Print the headings for Report #1, "Velocity Map for Inbound/Outbound Trajectories". 

10. Echo back the input values onto Report #1. 

11. Initialize the velocity map matrix coordinates to 1,1 (II = rows, or latitude increments; 
NN = columns, or longitude increments). Set the flag IPRINT to zero, which has the 
effect of causing the longitude increments to be printed as a subheading. Initialize the 
minimum AV hold variable (DVMIN) to 99999. 

12. Calculate the orbi tal velocity o f the Moon. 

YDLO =Y|xE/RREM 

13 Convert the initial map latitude and longitude from degrees to radians (ALAT and 
ALON). 

14. Initialize DVSIM at 99999. This variable will be used to hold the lowest AV for the 
trajectory, selected from all combinations of latitude and longitude. 

15. Initialize the current Total AV matrix cell (DELV) to 99999. 

16. Initialize with zeros the print line variables for the current cell on each report (PAGE1, 
PAGE2, PAGE3, PAGE4, PAGE5, PAGE6). 

17. Determine whether the current latitude for the map is less than the Lunar orbit inclina¬ 
tion. If the absolute value of the map latitude is greater than the inclination, then set the 
current latitude/longitude cell for all the report matrices to zero. Otherwise, perform the 
following iteration 41 times, incrementing the value of VELM by 100 feet per second 
each iteration, beginning at 800 ft/sec and ending at 4800 ft/sec. 

a. Call the in-program subroutine given the current latitude, longitude, and VELM 
to determine AV, flight times, and nodes. 

b. Leave the iteration if the Earth-to-SOI flight time is near zero. 

c. If the total AV for this iteration of VELM is the smallest so far for the current 
latitude/longitude, then store the following information in the current long¬ 
itude/latitude cell of the appropriate report matrix: 

- Total AV (DVTOTAL) 

- SOI AV (DELVEL) 

- Velocity at SOI for the SOI-to-Moon trajectory (VELMOUT) 

- Earth and Lunar ascending/descending nodes (ANODEE, ANODEM) 

- Earth and Lunar circular orbit AV’s (DVCIRE, DVCIRM) 

18. Store the matrix column heading (longitude increments) for use in all the output reports. 
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19. If all the matrix columns have been processed for this latitude, and if this has been the 
first row of the matrix, print the matrix column headings for Report #1. 

20. If all ten matrix columns have not been processed for this latitude, increment the column 
counter by one and increment the longitude by the amount specified by the user in the 
inputs. Return to step 15 to process the next longitude/latitude combination. 

21. If all ten matrix columns have been processed for this latitude, print the matrix row for 
this latitude for Report #1. 

22. If all 19 matrix rows have not been processed, increment the row counter by one, 
increment the latitude by the amount specified by the user in the inputs, reset the column 
counter to one, reset the longitude to the initial input longitude, and return to step 15 to 
process the next longitude/latitude combination. 

23. If all 19 matrix rows have been processed, continue with the following steps. 

24. Set VELM to the SOI velocity associated with the minimum total AV. 

25. Set the SOI latitude (ALAT) and longitude (ALON) to those associated with the 
minimum total AV (ALATMIN, ALONMIN). 

26. Call the in-program subroutine to determine the orbital characteristics associated with the 
minimum total AV. 

27. Print the remaining five reports. 
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Appendix E - In-Program Subroutine Description 
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This section is a line-by-line description of the lines of code in the subroutine that is imbedded in 
program PLANECHG, beginning at line #25. Refer to Appendix B for the actual code listing, 
and figures J-2 through J-4 for supporting illustrations. 


I. Establish the key Earth, Moon, and Sphere of Influence distances and angles. 

A. Calculate the distance from the Moon to the SOI point. 

1. Define variables for the sine and cosine of the SOI latitude and longitude 
(COSALAT, COSALON, SINALAT, SINALON). 

2. Given a point on the SOI defined by the latitude and longitude, define the 
angle $ between the Moon-to-Earth line and the Moon-to-SOI line. 

cos($) = cos(lat) * cos(lon) 

3. Define variables for the sine and cosine of the Lunar orbit inclination 
(COSAINC and SINAINC). 

4. Calculate the ratio of Earth-Moon distance to Moon-SOI distance 
(RMR = REM / RM). 

a. RRE 2 = RRM 2 + RREM 2 - 2 * RRM * RREM * cos(#) 

b. [lE/RRE 2 = (J.M/RRM 2 RRE 2 /RRM 2 = pE/pM 

(because |i/R 2 = gravitational acceleration, and by definition, at the 
sphere of influence, acceleration towards the Earth equals 
acceleration towards the Moon). 

c. RRE7RRM 2 =pE/|xM= 1+RREM7RRM 2 -2*(RREM/RRM)*cos(4>) 
= 1 + RMR 2 - 2*RMR*cos(«) 

d. 0 = -(|iE/|iM) + 1 + RMR 2 - 2*RMR*cos(4>) 

= RMR 2 - 2*RMR*cos($) +1 - (pE/^lM) 

e. For the quadratic formula: 

a = 1, b = -2*cos(§), c = 1 - (|iE/pM) 

RMR = (2*cos($) ±V4cos 2 - 4(1 - pE/(iM)’ ) / 2 
RMR = cos(§) +Ycos 2 -1 + liE/pM 1 

5. Determine the distance from the Moon to the SOI. 

RRM = RREM/RMR 
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B. Determine an Earth-based rectangular coordinate system such that the X-direction 
is on the line from the Earth to the Moon; the Y-direction is in the direction of the 
Moon’s orbit; and the Z-direction is perpendicular to X and Y in right-hand 
coordinates. Project the SOI point onto the X-Y plane (Earth-Moon plane) in 
order to calculate the X- and Y-coordinates. 

a. Determine the distance from the Moon to the X-intercept. 

XXM = RRM * cos($) = RRM * COSALAT * COSALON 

b. Determine the distance from the Earth to the X-intercept. 

XX = RREM - XXM 

c. Determine the Y-coordinate. Note that longitude increases in the negative 
Y direction. 

YY = -RRM * COSALAT * SINALON 

d. Determine the Z-coordinate. 

ZZ = RRM * SINALAT 

C. Calculate the distance from the Earth to the SOI point. 

RRE =1(XX 2 + YY 2 + 72? 

D. Identify the angle between the Earth-to-Moon line and the Earth-to-SOI line. 

a. cos(ANGA) = XX / RRE 

b. sin(ANGA) =lYY 2 + ZZ 2 ' / RRE 

c. ANGA = tan'(sin(ANGA) / cos(ANGA)) 

E. Determine the Earth-based latitude and longitude of the SOI point. 

a. ALONX = tan' 1 (YY / XX) 

b. ALATX = tan' 1 (ZZ / YxX 2 + YY 2 ') 

F. Determine the Earth ascending node for the orbit in which the trans-SOI bum will 
occur. 


a. Given the latitude of the bum at Earth (negative SOI latitude) and the 
inclination of the Earth orbit, spherical coordinate trigonometry states that 
the longitude from the node to the bum point is 

x> = sin' (tan(-ALATX)/tan(AINCE)) 

b. Longitude of the bum point (measured from the Earth-Moon line) is 
ALONX. 
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c. The node is calculated by subtracting the longitude in (a) above from the 
longitude in (b). 

ANODEE = ALONX - sin’ 1 (tan(-ALATX) / tan(AINCE)) 

II. For the SOI-to-Moon trajectory, calculate velocity at the SOI point, AV at Lunar circular 

orbit, Lunar ascending/descending node, and time of flight. 

A. Calculate flight path angle, AV at Lunar circular orbit, and time of flight. 

a. Call the subroutine GAMACALC to determine the cosine of the flight 
path angle, Lunar perigee velocity. Lunar circular orbit velocity, and time 
of flight. 

b. The cosine of the flight path angle will be between zero and one, 

corresponding to angles between zero and 90. In reality, Moon-to-SOI 
flight path angles will range from zero to 90, but SOI-to-Moon angles will 
range from zero to -90. Gamma must be derived from cos(GAM) and 
multiplied by -MOOD (from the input: +1 for outbound and -1 for return) 
to get the correct flight path angle. The arctan function is used instead of 
arccos to evaluate gamma because of that function’s better debugging 
diagnostics capabilities. r—_7 

GAMAM = -MOOD * tan'(il - COSGAM 1 / COSGAM) 

B. Identify the angle between the Moon-to-Earth line and the Moon-to-SOI line. 

a. Calculate $ (1) for node option 1. 

$ (1) = sin '(sin(latitude) / sin(inclination) 

$ (1) = sin' 1 ((ZZ / RRM) / sin(AINCM)) 

b. Calculate # (2) for node option 2. 

<t> (2) = II - f (1) if ZZ is positive, or 

$> (2) = -$(1) - II if ZZ is negative, or 

$ (2) = II if ZZ is zero. 

C. Determine the azimuth of the SOI point. 

AZMM = tan' (cot(AINCM)/cos(3>(NP))) 

= tan' 1 (cos(AINCM) / (sin(AINCM) * cos(<? (NP))) 
where NP is the node option from the input. 

D. Determine the Lunar ascending node. 

a. Calculate v (upsilon) for a Lunar orbit that has inclination AINCM and 
that passes through the SOI point. 

AVM = tan 1 (sin(ALAT) / (cos(*) * tan(AINCM))) 

= tan' 1 ((tan(ALAT) * cos(ALAT)) / (cos(<f) * tan(AINCM))) 

= tan' 1 ((tan(ALAT) / tan(AINCM)) / (cos(£) / cos(ALAT))) 
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b. Determine the node. 

ANODEM = ALON - AVM 

E. Determine the components of the velocity vector at the SOI point. 

a. Call the subroutine VELTRANS to determine the X-, Y-, and Z- 
components of the velocity vector, from the viewpoint of the moon. 

b. Add the X- and Y-components of the Moon’s velocity to determine the 
SOI velocity from the viewpoint of the Earth. 

VXM = VXM + XDLO VYM = VYM + YDLO 

ITT For the Earth-to-SOI trajectory, calculate velocity at the SOI point, AV at Earth circular 
orbit, and time of flight. 

A. Modify the velocity at SOI, if necessary, to bring the trajectory up to at least the 
minimum velocity required to intercept the SOI. 

a. Calculate the vis-viva parameter QQEMIN for the Earth-to-SOI trajectory. 

QQEMIN = 2 - (RRE / a) where a = (RPE + RRE) / 2 
QQEMIN = 2 - ((2 * RRE) / (RPE + RRE)) 

QQEMIN = 2 * (1 - (RRE / (RPE + RRE))) 

QQEMIN = 2 * (RPE / (RPE + RRE)) 

b. Calculate the minimum velocity required, such that apogee is just at SOI. 

QQEMIN s (V 2 * RRE) / CMUE _ 

.-. VELEMIN =Y(QQEMIN * CMUE) / RRE 1 

c. Increase the velocity slighdy so that the calculations converge. 

d. Compare the most recently calculated SOI velocity (VELE) with the 
minimum velocity. If it is too low, bring it up to the minimum. 

B. Temporarily store (until step IV-D) the calculated Earth-to-SOI time of flight. 

C. Further increase the velocity at SOI, if necessary, to meet the Earth-to-SOI time 
of flight requirement (but capping the AV at 500 ft/sec). Calculate the new SOI 
flight path angle, AV at Earth circular orbit, and time of flight. 

a. Call the subroutine GAMACALC to determine the time of flight from 
Earth perigee to SOI (TIMEE1), given the velocity VELE. 

b. Increase the velocity by 10 ft/sec, and call GAMACALC again to 
determine the new time of flight (TIMEE2). 
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c. For the Earth-to-SOI leg, determine the ratio of time-of-flight shortfall to 
the time-of-flight increase that was due to the 10 ft/sec increase in 
velocity. 

Shortfall/Increase = (FTIME-TIMEE1 -TIMEM) / (TIMEE2 - 
TIMEE1) 

Apply this ratio to the increased velocity of 10 ft/sec to estimate the AV 
necessary to meet the required time of flight. 

DELVELE = 10 * (FTIME - TIMEE1 - TIMEM) / (TIMEE2 - 
TTMEE1) 

d. Cap this AV at 500 ft/sec, and add it to the original velocity. 

e. Call the subroutine GAMACALC to determine the new Earth-to-SOI time 
of flight for the revised velocity. 

D. Compare the total Earth-to-Moon time of flight to the required time of flight. If 

the total is short, return to step HI-C. 

E. Determine the components of the velocity vector at the SOI point. 

a. Calculate the flight path angle at SOI (GAMAE) for the Earth-to-SOI 
trajectory (see step II-A-b for a discussion of this). 

b. Determine the longitude of the SOI from the Earth’s viewpoint. 

ALONE = ALONX +180 (since zero longitude is on the Earth- 
Moon line on the side of the Earth away from the Moon). 

c. Determine the latitude of the SOI point. 

ALATE = tan 1 (ZZ / 7XX 2 + YY 2 ) 

d. Calculate the azimuth of the SOI point 

AZME = tan' (cot(AINCE) / cos(PHIE)) 

AZME = tan' (1 /(tan(AINCE) * cos(PHIE))) 

where PHIE = II - sin'((ZZ/RRE)/sin(AINCE)) - see step H-B-b 

for a discussion of this. 

e. Call the subroutine VELTRANS to determine the X-, Y-, and Z-compo- 
nents of velocity at SOI for the Earth-to-SOI trajectory. 

IV. Combine the SOI-to-Moon calculations with the Earth-to-SOI calculations to identify 
total flight characteristics. 
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A. Calculate total AV 

a Calculate AV at the SOI point required to patch the Earth-to-SOI 
trajectory into the SOI-to-Moon trajectory. (This value is calculated by 
determining the difference between each of the components of the two 

"dELVEL =Y(VXE - VXM) 2 + (VYE - VYM) J + (VZE - VZM) J ' 

b. Calculate AV at Lunar perigee. 

DVCIRM = VPM - VCIRMP 


Calculate AV at Earth perigee. 

DVCIRE = VPE - VC IRE 


d Save the old total AV. 

DVTSAV = DVTOTAL 


B. 


Calculate the new total AV. 

DVTOTAL = DELVEL + DVCIRM + DVCIRE 


Determine the distance from the Moon to the Earth for this trajectory just 
calculated. Use this value in subsequent iterations of the in-program subroutine, 

if necessary. 


Determine the time of arrival at SOI from Earth, in centuries since the 
year 2000 AD. 

TIM = Departure date + Time of flight - 2000 


where Departure date = TIMJ / 36525 days per century 
Time of flight = TIEM hrs / 876600 hrs per century 
The year 2000 = Day 2451545 / 36525 days per century. 


Call the subroutine POSVELMO to determine the Moon’s distance from 
the Earth and velocity at the time of arrival at SOI from the Earth. 


c. 


Convert Lunar distance from Earth radii to feet. 
RREM = RREMER * 20295741 ft/radii. 


C. Compare the value of Earth-to-SOI time of flight saved in step DI-B above with 
the new value of time of flight calculated in step m-C above. If the difference 
has not yet converged (to less than 0.1 hour), reiterate this in-program subroutine 
beginning at step I-A. Otherwise, continue. 
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If the newly calculated total AV is the smallest AV calculated so far then store it 

Liable DWM ed A S , 01 AV * the SmaUeSt CaIculated so far > then store it in 
the variable DVSIM. Also store its corresponding SOI longitude (ALONSIM) 

fvELMSIM) ALATSIM) ’ md S01 velocit y for ‘he SOI-to-Moon trajectory^ 
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Appendix F - Description of Subroutine GAMACALC 
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1 . 


The subroutine GAMACALC receives the parameters orbital perigee radius tRPX'* nrhit^i 

^ingmbhed^l^X) 01 !^ al Ia f r ** (RRX) ’ 30(1 the gravitational parameter of the body 

path a^Se rCOST AMX^ 1 retUms tlme of flight H1MBX). the cosine of the flight 

path ^gle (COSGAMX), velocity at penapses (VPX), circular velocity at periapses (VCIRX) 

an mdicator descnbmg whether the orbit is elliptical or hyperbolic (TRAJ). 

Initialize indicators to presume an elliptical orbit. (IHYPER, TRAJ). 

Calculate the vis-viva parameter 

QQX = (RRX * W 2 ) / CMUX. 

If the orbit is just barely hyperbolic (QQX is within one-millionth of 2), force the 
calculation to consider it elliptical (reduce QQX to just under 2). 

!rri C n° rbit “ ^ h3 y erbo u lic ’ reset the indicators to show this. Print a message on the 
screen announcing a hyperbolic orbit. * 

Calculate the semi-major axis of the orbit 
AAX = RRX / (2-QQX). 

If the semi-major axis is very large (greater than 10'*) or very small (less than -10") the 
orbit is trapped near a parabolic trajectory. Make it hyperbolic: 

AAX ~ -10'\ 


3. 


4. 


6 . 


7. Calculate the orbit eccentricity 

EEX = 1 - (RPX/AAX). 

8. Calculate the semi-latus rectum 

PPX = AAX * (1 - EEX 2 ). 

9. Calculate the flight path angle 

a. The angular momentum of the orbit at perigee is 

HP = RPX * VELPERIGEE * cos(GAMAX). 

But at perigee, GAMAX is zero, so 
HP = RPX * VELPERIGEE. 

b. At SOI, angular momentum is 

HX = RRX * VELX * cos(GAMAX). 

c. Angular momentum is constant along a given orbit so 

HP = HX 

RPX * VELPERIGEE = RRX * VELX * cos(GAMAX) 
GAMAX = cos '((RPX * VELPERIGEE) / (RRX * VELX)) 
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d. 


10 . 

11 . 

12 . 

13. 

14. 


15. 


f. 


™ c vel ^ E S e E e 4cMux * RAPOGEeT 7(AAX*RPX?. 

Since RAPOGEE / A AX = (1 + EEX). then_ 

VELPERIGEE =TCMU X * (1 +jpX) / 

RPX * VELPERIGEE =1RPX * (1 + EEX) * cMUX 1 . 

Substituting (d) into (c) abo ve yields - VELX'rt 

GAMAX = cos (t^RPX * (1 + EEX) * CMUX) / (RRX VELX)) 

Substituting from (2) above yields— 

GAMAX = cos' 1 tfRPX * (1 + EEX)) / (RRX QQX)). 


Calculate the p erigee velocity __ 

VPX =^CMUX * (1 - EEA) / RKL 

Calculate the circ ular velocity at perigee 
VCIRX =YCMUX / RPX. 

Calculate the true anomaly 

THETAX = cos* , (((PPX / RRX) -1) / EEX) 

Calculate the eccentric anomaly rrunTAY'i^^ 

AEX = cos '((EEX + cos(THETAX)) / (1 + EEX * cos(THETAX))). 

Compare the orbital radius at SOI (RRX) to the apogee radius (AAX * U + EEX)). If the 
orbital radius at SOI is greater than the apogee radius, print a message on e 
indicating the difference in nautical miles. Also display the following. 

a. SOI latitude (ALAT) 

b. SOI longitude (ALON) 

c. vis-viva parameter (QQX) 

d. semi-major axis (AAX) 

e. eccentricity (EEX) 

f. flight path angle (GAMAX) 

g. perigee velocity (VPX). 

Calculate the time of flight. 

a . If the or ^^ e ^^^3 /CMU ^ * (aex - EEX * sin(AEX)). 

b. If the or ^JJ^ h ^^^3 /( 3^*(EEX*sinh(EEX) - log(cosh(EEX) + 

sinh(EEX))). 
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Appendix G - Description of Subroutine POSVELMO 
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The subroutine POSVELMO receives the parameter TIME (number of Julian centuries from the 
year 2000) and returns the moon’s position and velocity at that time. Specifically, it returns the 
moon’s distance from the Earth’s center, in Earth radii (RRM); velocity in the x-direction (along 
the Earth-Moon line), in feet per second (XDLO); and velocity in the y-direction (direction of 
the moon’s orbit), in feet per second (YDLO). 

1. Calculate a fraction of TIME that represents half a second. 

DELT = 0.5 seconds / (36525 days/century * 24 hrs/day * 3600 seconds/hr) 

= 1.58440E-10 centuries. 

2. Call the subroutine MOON with the parameter (TIME minus DELTA) to determine the 
moon’s distance in Earth radii at half a second before TIME. This distance is RM1. 

3. Call the subroutine MOON with the parameter (TIME plus DELTA) to determine the 
moon’s distance in Earth radii at half a second after TIME. This distance is RM2. 

4. Calculate the velocity of the moon in the x- (radial) direction. 

XDLO = [(RM2 Earth radii - RM1 Earth radii) / (1 sec)] * 20,925,741 ft/radii. 

5. Calculate the average radius of the Lunar orbit during the one second centered on TIME. 

RRM = (RM1 + RM2) / 2. 

6. Determine the radius of the Lunar orbit from the Earth-Moon baricenter. 

RRMB = RRM - 0.7412789 Earth radii. 

7. Calculate the moon’s velocity in the direction of its orbit. 

a. Moon’s apogee (Apo) = 62.83308 Earth radii. 

b. Moon’s perigee (Per) = 55.68264 Earth radii. 

c. Eccentricity (e) = (Apo - Per) / (Apo + Per) = 0.06033. 

d. Semi-latus rectum (p) = Apo(l -e 2 ) = 62.60439 Earth radii 

= 1,310,038,967 feet. 

e. Earth’s gravitational parameter (mu) = 1.407646822E+16 ff/sec 2 . 

f. Angular momentum (h) =Ymu * p' 

= 4.29427E+12 ft 2 /sec * (1 earth radii / 20,925,672.57 ft) 

= 205,215.4 ft*Earth radii/second. 

g. Y-velocity (YDLO) = h / RRMB = 205,215.4/RRMB. 
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Appendix H - Description of Subroutine MOON 
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The subroutine MOON receives the parameter T (number of Julian centuries from the year 2000) 
and returns the approximate location of the moon in geocentric coordinates at that time. 
Specifically, it returns the right ascension of the moon (RAM), declination of the moon 
(DECM), and distance to the moon in Earth radii (RM). The formulae are from The Astronomi¬ 
cal Almanac For The Year 1984, page D46. All degrees are converted to radians with the 
conversion factor C = re/180. 

1. Calculate the ecliptic coordinates of the moon. 

a. Moon’s ecliptic longitude 

LAM = 218°.32 + 481267°.833T 

+ 6°.29 * sin(134°.9 + 477198°.85T) 

-1°.27 * sin(259°.2 - 413335°.38T) 

+ 0°.66 * sin(235°.7 + 890534°.23T) 

+ 0°.21 * sin(269°.9 + 954397°.70T) 

- 0°.19 * sin(357°.5 + 35999°.05T) 

- 0°.ll * sin(186°.6 + 966404°.05T) 


b. Moon’s ecliptic latitude 

BETA = 5°. 13 * sin( 93°.3 + 483202°.03T) 
+ 0°.28 * sin(228°.2 + 960400°.87T) 

- 0°.28 * sin(318°.3 + 6003°. 18T) 

- 0°.17 * sin(217°.6 - 407332°.20T) 


c. Horizontal parallax 

PIE = 0°.9508 

+ 0°.0518 * cos(134°.9 + 477198°.85T) 

+ 0°.0095 * cos(259°.2 - 413335°.38T) 

+ 0°.0078 * cos(235°.7 + 890534°.23T) 

+ 0°.0028 * cos(269°.9 + 954397°.70T) 

d. Semi-diameter of moon’s orbit 

SD = 0.2725 * PIE 

e. Distance to the moon in Earth radii 

RM = 1 / sin(PIE) 

2. Form the geocentric direction cosines to rotate into geocentric coordinates. 

a. 1 = cos(BETA)cos(LAM) 

b. m = 0.9175*cos(BETA)sin(LAM) - 0.3978*sin(BETA) 

c. n = 0.3978*cos(BETA)sin(LAM) + 0.9175*sin(BETA) 

where 1 = cos(DECM)cos(RAM), m = cos(DECM)sin(RAM), n = SIN(DECM). 


3. Then: 

a. RAM = arctan(m/l) [right ascension] 

b. DECM = arcsin(n) [declination] 



The errors will rarely exceed 0.2 Earth radii in distance (RM), 0°.3 in right ascension (RAM), 
and 0°.2 in declination (DECM). 
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Appendix I - Description of Subroutine VELTRANS 
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The subroutine VELTRANS converts an orbital velocity vector into rectangular coordinates (see 
figure J-5). Parameters received by the subroutine are velocity (VEL), flight path angle 
(GAMA), azimuth (AZM), latitude above the Earth-Moon plane (ALAT), and longitude from 
the Earth-Moon line (ALON). A set of intermediate calculations is performed to express the 
velocity vector in terms of a radial component, a latitudinal component, and a longitudinal 
component. Each of these three components is further resolved into x-, y-, and z-components. 
Finally, all three x-components, all three y-components, and all three z-components are summed 
to provide the total x-, y-, and z-components of velocity. 

1. Conversion of velocity vector into spherical coordinates. 

From the geometry, the radial component of velocity, R, is calculated to be 
VEL * sin(GAMA). (See figure J-6). 

The component along the orbital path, r|, is 
VEL * cos(GAMA). 

This orbital path component of velocity is further resolved into a latitude component, 

LAT, and a longitude component, LON (see figure J-7). Again, from the geometry, 

LON = R$ * sin(AZM) and 
LAT = R* * cos(AZM). 

2. Conversion of spherical coordinates into rectangular coordinates. 

a. Conversion of radial component into rectangular coordinates. 

Refer to figure J-8. The projection of R onto the x-y plane is 
R * cos(LAT). 

The negative x-component of this is 
R * cos(LAT) t * cos(LON) 
so the x-component, X, is 

-R * cos(LAT) * cos(LpN). 

The negative y-component of R is 
R * cos(LAT) 4 * sin(LON) 
so the y-component, Y, is 

-R * cos(LAT) * sin(LON). 

The z-component, Z, is 
R * sin(LAT). 

b. Conversion of latitude component into rectangular coordinates. 

Refer to figures J-9 and J-10. LAT is perpendicular to the radial vector, R. A 
line in the z-direction that meets the tip of LAT and intersects the radius vector 
produces the angles a and b, where 
b = 90 - LAT and 
a + b + 90 = 180. Therefore, 
a = LAT. 

From the geometry, the z-component of LAT, 7.1-A T, is 
LAT * cos(LAT). 
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The projection of LAT onto the x-y plane is 
lAt * sin(LAT). (see # figure J-l 1). 

The x-component of this, XLAT, is 
LAT * sin(LAT) * cos # (LON). 

The y-component of this, YLAT, is 
LAT * sin(LAT) * sin(LON). 

c. Conversion of longitude component into rectangular coordinates. 

Refer to figures J-12 and J-13. LON is always parallel to the x-y plane, so the z- 
component of l6n, ZL6n, is always zero. Using the same proof as in (b) above, 
it can be seen that the angle between LON and the ^-component of LON is equal 
to LON. From the geometry, the x-component of LON, XLON, is 
LON * sin(LON). , 

The negative y-component of LON is 
LON * cos(LON|, 
so the y-component, YLON, is 
-LON * cos(LON). 

3. Sum of the rectangular coordinates. 

All of the x-, y-, and z-components are summed to provide the complete rectangular 
coordinates of the velocity vector. 

VXX = X + XLXT + XLON 
VYX = Y + YLXT + YLjbN 
VZX = Z + ZLAT + ZLON. 
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Appendix J - Illustrations 
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Figure J-1 











Figure J-7 
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Figure J-9 
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Figure J-10 


96 





9 : 




-X 


Figure J-13 
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